rm(list=ls())

# Load packages
library(foreign)
library(lmtest)
library(sandwich)
library(stargazer)
library(msm)
library(readxl)
library(ggplot2)
library(gridExtra)
library(readstata13)
library(haven)
library(dplyr)
library(forcats)
library(ggeasy)

# Load data
cep_all <- read_dta("Data/cep_all.dta")

# Define dosage range samples
dosage_ranges <- c(1, 2, 3, 4, 5)
sample_list <- lapply(dosage_ranges, function(i) subset(cep_all, treatment_dosage_coup >= -i & treatment_dosage_coup <= i))
names(sample_list) <- paste0("sample", dosage_ranges)

# Clustered SE function
cluster_se <- function(model, data, type = "HC1") {
  idx <- as.integer(rownames(model.frame(model)))
  vc  <- sandwich::vcovCL(model, cluster = data$cluster[idx], type = type)
  lmtest::coeftest(model, vcov. = vc)
}

# Reusable model block generator
run_models <- function(samples, formula, dv_name, table_title, dep_caption, column_labels, cov_label) {
  models_lm <- lapply(samples, function(s) lm(formula, data = s))
  models_cl <- Map(cluster_se, models_lm, samples)
  sample_sizes <- sapply(models_lm, nobs)
  r_squared <- sapply(models_lm, function(m) summary(m)$r.squared)
  mean_dv <- sapply(models_lm, function(m) mean(m$model[[dv_name]], na.rm = TRUE))
  
  stargazer(models_cl,
            type = "latex",
            title = table_title,
            keep = "^treat1$",
            column.labels = column_labels,
            covariate.labels = cov_label,
            dep.var.caption = dep_caption,
            dep.var.labels = dv_name,
            dep.var.labels.include = FALSE,
            add.lines = list(
              c("Obs", sample_sizes),
              c("R^2", round(r_squared, 3)),
              c("Mean DV", round(mean_dv, 3))
            ),
            align = TRUE)
}

# ------------------
# Table B1 - Education: College (5 & 6)
samples_college <- lapply(sample_list, function(s) subset(s, ed_levels %in% c(5, 6)))
tabb1_morehs <- run_models(samples_college,
           none_pol ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1+ as.factor(encuesta),
           dv_name = "none_pol",
           table_title = "Table B1: RD Estimates 1973 Coup on No Ideology by Education (College)",
           dep_caption = "Dependent Variable: No Ideology",
           column_labels = c("54 and 56", "53-57", "52-58", "51-59", "50-60"),
           cov_label = "Coup")

# Table B1 - Education: High School (4)
samples_hs <- lapply(sample_list, function(s) subset(s, ed_levels == 4))
tabb1_hs <- run_models(samples_hs,
           none_pol ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "none_pol",
           table_title = "Table B1: RD Estimates 1973 Coup on No Ideology by Education (High School)",
           dep_caption = "Dependent Variable: No Ideology",
           column_labels = c("54 and 56", "53-57", "52-58", "51-59", "50-60"),
           cov_label = "Coup")

print('Table B1 complete')
# ------------------
# Table B2 - Age < 55
samples_age_lt55 <- lapply(sample_list, function(s) subset(s, edad55 == 0))
tabb2_agelt55 <- run_models(samples_age_lt55,
           left ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "left",
           table_title = "Table B2: RD Estimates 1973 Coup on Left-wing Identification (Age < 55)",
           dep_caption = "Dependent Variable: Left Identification",
           column_labels = c("54 and 56", "53-57", "52-58", "51-59", "50-60"),
           cov_label = "Coup")

# Table B2 - Age ≥ 55
samples_age_gte55 <- lapply(sample_list, function(s) subset(s, edad55 == 1))
tabb2_agegt55 <- run_models(samples_age_gte55,
           left ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "left",
           table_title = "Table B2: RD Estimates 1973 Coup on Left-wing Identification (Age ≥ 55)",
           dep_caption = "Dependent Variable: Left Identification",
           column_labels = c("54 and 56", "53-57", "52-58", "51-59", "50-60"),
           cov_label = "Coup")

print('Table B2 complete')
# ------------------
# Table B3 - By Region
regions <- list(
  RM = subset(cep_all, region_rm == 1),
  South = subset(cep_all, region_south == 1),
  North = subset(cep_all, region_north == 1)
)

for (region_name in names(regions)) {
  samples_region <- lapply(dosage_ranges, function(i) subset(regions[[region_name]], treatment_dosage_coup >= -i & treatment_dosage_coup <= i))
  run_models(samples_region,
             left ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
             dv_name = "left",
             table_title = paste("Table B3:", region_name, "Region - RD Estimates 1973 Coup on Left-wing Identification"),
             dep_caption = "Dependent Variable: Left Identification",
             column_labels = c("54 and 56", "53-57", "52-58", "51-59", "50-60"),
             cov_label = "Coup")
}

print('Table B3 complete')
# ------------------
# Table B4 - By Period
samples_pre2014 <- lapply(sample_list, function(s) subset(s, encuesta_a >= 1994 & encuesta_a <= 2014))
samples_post2014 <- lapply(sample_list, function(s) subset(s, encuesta_a > 2014))

tabb4_pre2014<- run_models(samples_pre2014,
           left ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "left",
           table_title = "Table B4: RD Estimates 1973 Coup on Left-wing Identification (1994–2014)",
           dep_caption = "Dependent Variable: Left Identification",
           column_labels = c("54 and 56", "53-57", "52-58", "51-59", "50-60"),
           cov_label = "Coup")

tabb4_post204 <- run_models(samples_post2014,
           left ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "left",
           table_title = "Table B4: RD Estimates 1973 Coup on Left-wing Identification (Post-2014)",
           dep_caption = "Dependent Variable: Left Identification",
           column_labels = c("54 and 56", "53-57", "52-58", "51-59", "50-60"),
           cov_label = "Coup")

print('Table B4 complete')

